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ABSTRACT 

Finite-difference simulations of fluid dynamics and magnetohydro dynamics generally 
require an explicit diffusion operator, either to maintain stability by attenuating grid- 
scale structure, or to implement physical diffusivities such as viscosity or resistivity. 
If the goal is stability only, the diffusion must act at the grid scale, but should affect 
structure at larger scales as little as possible. For physical diffusivities the diffusion 
scale depends on the problem, and diffusion may act at larger scales as well. Diffusivity 
undesirably limits the computational timestep in both cases. We construct tuned finite- 
difference diffusion operators that minimally limit the timestep while acting as desired 
near the diffusion scale. Such operators reach peak values at the diffusion scale rather 
than at the grid scale, but behave as standard operators at larger scales. We focus on the 
specific applications of hyper diffusivity for numerical stabilization, and high Schmidt 
and high Prandtl number simulations where the diffusion scale greatly exceeds the grid 
scale. 



1. Introduction 

Fluid dynamics simulations usually use explicit diffusion operators, either to maintain stability 
or to model physical effects such as viscosity, resistivity, conductivity, or the diffusion of passive 
scalars. Virtually all astrophysical gas dynamics and MHD simulations rely on such diffusion 
operators for stability, physical effects, or both. We here consider how to design such operators 
such that they have the desired behavior at the diffusion scale and larger scales, while still restricting 
the numerical timestep as little as possible. The timestep depends inversely on the strength of the 
diffusion 

At = Ax 2 /2v. (1) 

Classical diffusion operators such as Laplacian viscosity (z^V 2 ), or fourth or sixth-order hyperdif- 
fusivities (z^V 4 or fgV 6 ), reach their maximum values at the grid scale, but act at the larger 
diffusion scale where the effective diffusivity is lower. The operators designed here reach their max- 
imum value at the diffusion scale rather than at the grid scale, so that they limit the timestep no 
more than necessary. 
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In a spectral code, the diffusive terms are linear and can thus be handled spectrally without 
limitation on the timestep. For example, let a field evolve as c^V = A — z^V 2 V, where A denotes the 
non-diffusive terms. In Fourier space, dfV = A — vk 2 \ The solution, with A constant throughout 
the interval At, is 



V(At) 



V(0) + 



vk 2 At 



I) 



-vk 2 At 



(2) 



When evolved in Fourier space, the diffusivity operator is stable for any value of vk 2 At, whereas 
in physical space, instability occurs if vk 2 At > 2 (restating Eq. [I] in terms of wavenumber). 

However, finite difference codes do have advantages that make them worth pursuing: they 
use fewer floating point operations per grid point; they can be more easily parallelized without 
the all-to-all communications required for Fourier transforms; they are not restricted to periodic 
boundary conditions; and they handle discontinuous jumps more robustly. 

The Navier-Stokes equation can include a number of different types of diffusion operators: 



d t V = V- VV-p _1 VP + i/ 2 V 2 V-i/ 4 V 4 V + i/ 6 V 6 V 

- v' 4 (di + d 4 y + d 4 z )V + u'M + d 6 y + 6 Z )V - v D D[V) 



(3) 



where the term is the usual Laplacian physical viscosity, the v n and v' n terms are nth-order 
hyperviscosities, the term udD(V) is a customized diffusion operator. 

Either the sixth-order hyperdiffusivity term or t he physical diffusivit y can maintain numerical 
stability. The hyperdiffusivity has been advocated ([Brandenburg] l2003l ) because it preferentially 
diminishes the high-wavenumber structure without modifying low- wavenumber structure. If the 
problem does require true physical diffusivities, we still want to consider use of a customized 
operator. This would reduce excess diffusion at scales well below the diffusion scale. Such excess 
diffusion limits the timestep without further modifying the solution as no structure exists at those 
scales. 

To date the focus in the study of extensions to nume rical diffusion has been on such hyper- 
diffusivities (e.g. iBorue & Orszaglll995l : lBrandenburg)l2003l ). as exemplified by the hyperviscosities 
described in equation ([3]). However, the degrees of freedom available in the finite difference coeffi- 
cients can be used instead for different goals. 

In this paper we describe methods for customizing diffusion operators that can be used to design 
operators that protect the timestep while either minimizing diffusion or reproducing the physical 
diffusion operator at low wavenumber as well as possible. These methods can also be used to design 
di fferent diffusion ope rators for other purposes. These methods rely on the tuning techniques used 
by lMaron et al.l (|2008l ) for improving the high-wavenumber accuracy of finite-difference derivatives. 



In §[2] we summarize the constraints that lead to the need for tuning finite difference operators. 
We then describe operators suitable for implementing both numerical (§ [3|) and physical (§ 2]) 
diffusivities, and we summarize our results in § [5j 
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2. Tuning Finite Difference Operators 

Let us consider the question of how to tune a general, symmetric, finite-difference operator, 
since all diffusion operators must be symmetric. We follow the tr e atmen t of the tuning of anti- 



symmetric operators such as first derivatives given in iMaron et al.l (|2008l ) . Such tuning allows us 
to customize the wavenumber spectrum of the operator to meet the needs of the problem at hand, 
rather than relying on simple analytic forms. This allows us, for example, in a problem with large 
Laplacian diffusivity, to trade small deviations from Laplacian behavior at low wavenumber for large 
gains in the timestep, by limiting the diffusion at high wavenumber. The deviations from Laplacian 
at low wavenumber can be maintained at levels small enough to not affect realistic simulations. We 
use both analytic solutions and numerical optimization to improve the spectral performance of the 
operators. 

As an example of finite difference representations of symmetric operators we examine second 
and fourth derivatives. Define a function fj(xj) on a set of grid points Xj = j, with j an integer. 
Then construct a finite difference operator for the second derivative /t 2 ' by sampling a stencil of 
grid points with radius S. Without loss of generality, we center the operator on j = and use a 
grid interval of Ax = 1. The familiar result for a second derivative on a radius- 1 stencil is 

0*/(x)U =o ~-2/o + (4) 

which is obtained from fitting a polynomial of degree 2 to fj. For a fourth derivative, we can fit a 
degree 4 polynomial on a radius-2 stencil, 

d 4 J(x)\ x=0 ~ -6/o + 4(/_i + /i) - (f- 2 + /a). (5) 

In general, a symmetric operator on a stencil of order S can be represented as 

s 

i=i 



Consider the value of the finite-difference operator at x = for a Fourier mode / = cos(7rA:x). 
(Sine modes can be ignored because they don't contribute to the second derivative at x = 0.) 
The wavenumber k is scaled to grid units so that k = 1 corresponds to the maximum (Nyquist) 
wavenumber ir(/S.x)~ l expressible on the grid. The analytic value for the second derivative is —ir 2 k 2 , 
whereas the finite difference operator (eq. [6]) gives 

s 

f [2] ~ m + 2^ rrij cos(vr jk) = -D(k) (7) 
j=i 

This defines a function D(k) that, when positive, acts as a diffusion applied to f(x), because the 
Fourier modes of / evolve as d t f = —V£,D{k)f^ where V£> is a viscosity-like parameter that sets the 
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level of diffusion. The maximum diffusive timestep is given by the inverse of the maximum value 
of D(k) over < k < 1. 

1 

< u D max[D(k)Y (8) 

Ideally, D(k) should scale as (irk) 2 for k < k^ and should be constant for k > k^. The focus of this 
work is on customizing the form of D{k) so as to to increase the maximum diffusive timestep, and, 
in the case of hyper diffusion, also minimize low-A; diffusion. 

Figure Q] shows D(k) for finite-difference stencils of radius S = 1 (second order) and 5 = 3 
(sixth order), in comparison to the analytic value, demonstrating how higher order more closely 
mimics the analytic function. The coefficients of these functions are listed in Table [TJ The operator 
D{k) can be Taylor expanded in the form 

D{k) = D + D 2 k 2 + + D 6 k 6 .... (9) 

An operator that reproduces d 2 for all k would have D 2 = n 2 and all the rest of the coefficients 
D n = for n ^ 2. The radius- 1 stencil (Eq. 2|) has Do = and D 2 = ir 2 , but the higher order 
coefficients are unconstrained, while the radius-2 stencil (Eq. H|) sets Do = D± = and D 2 = ir 2 . 

The usual way to evaluate an operator for a function such as d 2 is to fit a maximal order 
polynomial to the points in the stencil. This yields equations that can be inverted to find the rrij 
coefficients in Equation [7] for stencil radius S 

s 

D = -m -2^m 9 (10) 

9=1 

D P = -^(-D P/2 I> P "V 

9=1 

for even p. We call such operators polynomial-based operators. The form of these operators is shown 
in Figure [2j We may, however, use the available degrees of freedom in different ways. The goal 
of fitting a high-order polynomial to the derivative is to have high accuracy at high wavenumber. 
However, that may actually contradict the goal of protecting the timestep while either minimizing 
diffusion or reproducing the physical diffusion operator at low wavenumber. Instead, we can use 
the available degrees of freedom to directly address these requirements. 

As an example, for a Kolmogorov cascade, the diffusive scale \ u and the viscosity v scale as 
\ v ~ i/ 3 / 4 . The viscosity can be made sufficiently large that \ v is substantially larger than the grid 
scale, and the velocity profile will be smooth at smaller scales. The cascading energy is elminated 
at the diffusive scale, so there is no need for higher diffusivity at smaller scales, yet because of 
the form of the Laplacian diffusivity operator, the diffusivity increases all the way down to the 
grid scale. This excess diffusivity is unnecessary, and in fact is a liability because it restricts the 
timestep. 
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Laplacian diffusion operators, or steeper operators such as hyperdiffusivities, rise in amplitude 
monotonically all the way to the Nyquist wavenumber k = 1. A large value of -D(l) requires a small 
timestep, but is unnecessary because energy cascading from higher scales is removed earlier at the 
diffusion scale kd- D{k) need only have enough presence above kd to diffuse any Fourier modes that 
might arise there. For k > kd, it can be as large as it is at kd with no additional timestep penalty 

This allows us to specify a strategically chosen diffusion operator that satisfies these require- 
ments. Such an operator should rise through the diffusion range kd, but then flatten out and 
merely remain positive at k > kd- Since this constraint is much less critical than having controlled 
diffusion at low k, the low-A; range of D(k) should receive a higher priority in the optimization than 
the high-A; range. 

The behavior of the diffusion function at kd critically determines its effect on the solution. 
The diffusion must act above both the Nyquist scale k^Y = lj for stability, and at the resolution 
scale, to damp modes with wavenumbers too large to be accurately captured by the finite difference 
scheme. The resolution scale depends on the details of the method. H owever, the common choice 



of radius-3 stencils have a resolution scale k = 1/2 (jMaron et al.ll2008l ). so we choose in this work 



to use a diffusion scale kd = 1/2, and examine diffusion functions normalized to D{kd) = 1. 

We note in passing that if the diffusivity is weak enough to not limit the timestep, it does not 
have to be applied every timestep. The diffusivity can instead be applied onc e every N timestep s 



with a value of of v that is N times as large, for reasonable values of N (jMaron et al.l 120081 ) . 
However, the enhanced diffusivity may then be large enough that the flat diffusivities described 
in this paper are required to protect the timestep. This yields a computational savings from not 
having to calculate the diffusion operator every timestep. 

In the next two sections we describe how we perform the tuning and give some useful examples 
of operators for both hyperdiffusivity and Laplacian diffusivity. 



3. Timestep-friendly hyperdiffusion 

In situations where one wishes to maximize the scale range, and where the diffusion-scale 
dynamics don't affect larger scales, one can fruitfully use a diffusion operator that rises more 
rapidly with k than a Laplacian. Writing the diffusive terms from equation [3] in Fourier space, 

8 t Y = -v 2 k 2 V - u^V - v 6 k 6 V - v' A (k A x + k* + kj)V - v's(kl + k & y + k 6 z )V - v d D{k)\ (11) 

one sees that hyper diffusive terms such as those proportional to 1/4 and uq have a steeper dependence 
on the wavenumber k than the Laplacian diffusivity proportional to v<i . 

However, energy cascading from larger scales dissipates at the diffusion scale kd < 1, so in- 
creasing the diffusion at k > kd further is unnecessary. A customized diffusion operator D(k) can 
be constrained to have a similarly steep k dependence at low k, but to then flatten at the diffusion 
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scale kd, so that the value at the Nyquist wavenumber D(l) is not markedly higher. Since the 
timestep is limited by the maximum diffusivity on the grid at any scale, limiting the value of D at 
small wavenumber protects the timestep. 

We note in passing that the ^4 term in Equation (jllf) contains two successive Laplacians and 
therefore two rounds of finite differences, whereas terms such as {&£ + dy + df) and D (k) involve 
only one round of finite differences, and are therefore favored for their execution speed. Also, 
the diffusion function for V 4 has a greater value in the high-/c "corners" of Fourier space than 
(<9 4 + dy + <9 4 ), and hence a smaller maximum diffusive timestep, and so for this reason as well, 
operators such as V 4 and V 6 are disfavored. 

We begin by considering diffusion operators wit h a stencil radius S = 3, such as are used in 



the hyperdiffusion implemented in the Pencil code (IBrandenburg fe Doblerl I2002I ). In this case, 
one can analytically construct a one-parameter family of functions parameterized by the degree of 
diffusivity D{\) at k = 1. As before, we take -D(O) = 0, the diffusion scale = 1/2, and normalize 
D(k) so that .0(1/2) = 1. Inverting Equation [7] with these conditions yields 

m = ~ - \d{1) (12) 

mi = \ + ~h D{l) (13) 

™2 = \ ~ \d{1) (14) 
^3 = -l + ^D(l) (15) 
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(16) 



W(i-^). (17) 



This also implies that 

D.-^( 1 — 1 - 

16 

The magnitude of D4 is inversely related to the sharpness of the hyper diffusive filter. Sharpness of 
hyperdiffusivity is usually measured by giving the index of the scaling with wavenumber k at low k. 
In this example, however, all the functions we present have k scaling at low k and are normalized 
at kd, so the lower the value of the coefficient -D4 of the k term, the sharper the hyperdiffusivity. 

The free parameter D(l) traces the maximum diffusivity, at least in the regime D{1) > 2, as 
shown in Figure [3l and thus determines the timestep. This can be demonstrated by differentiating 
D(k) (Eqs. [7] and [2]) and showing that in this regime, D'(k) > 0, so D{k) monotonically increases 
between < k < 1. For 1 < D(l) < 2, the maximum diffusivity is not much greater than the 
diffusivity at k = 1. For example, for D(l) = 1.5, the maximum diffusivity is D = 1.63 at k = 0.762. 
The useful range for D(l) is 1 < D(l) < 8 because the case D(l) = 8 corresponds to the operator for 
<9 6 , which represents the S = 3 hyperdiffusivity operator that is least diffusive at low wavenumber. 
Choosing D(l) > 8 results in D(k) < for some value of k < kd- The choice D(l) = 8 corresponds 
to the standard hyperdiffusivity used in the Pencil code. We further find that the best choice to 
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minimize diffusivity at low k is given by requiring that D"(0) = 0. If D"(0) < 0, then D{k) < at 
low k, and hence unstable, while for D"(0) > 0, it is more diffusive at low k than it could be. 

The goals of minimizing the \ow-k diffusivity and protecting the timestep at high k are at odds 
if one normalizes the diffusion magnitude at the diffusion scale D{kd) = 1. The tradeoffs can be 
seen by considering the behavior of D(k) as D{1) is increased (Fig. [3]). Reducing the diffusion for 
k < kd requires increasing the diffusion for k > kj, and vice versa. Although the values are much 
larger at high wavenumber, the percentage changes are actually similar in the two regimes (see 
Tab. [2] for the low wavenumber values). 

Note that the function with D(l) = 4 is the standard hyper diffusion with stencil radius 5 = 2, 
which is the most hyper diffusive 5 = 2 operator. Adding one extra free parameter by moving to 
stencil size 5 = 3 allows both low and high wavenumber diffusivity to be tuned, but we cannot 
decrease both simultaneously. However, adding another free parameter by using stencil size 5 = 4 
does allow both to be decreased (Fig. [2]). As an additional example, note that the dotted line in 
Fig. [3] shows the 5 = 1 diffusivity, while the 5 = 3 result with D(l) = 1.5 has decreased diffusivity 
for both low and high k. Simultaneously decreasing the diffusion for k < kd and k > kd while 
maintaining a constant diffusion at k = kd thus clearly requires more than one free parameter. 
With two or more free parameters, we can simultaneously satisfy both goals. 

Extending the stencil size to 5 > 3 allows us to further optimize the diffusion function. We 
now have multiple ways in which we could proceed. We choose to use numerical optimization to 
derive tuned 5 > 3 operators based on the following conditions: normalize the diffusion spectrum 
to D(kj) = 1; insist that < D (k) < 5 for some chosen value of the constant 5; set D (0) = 0; and 
require that D monotonically increase for k < k^- Within these constraints, we maximize D'(kd), 
which measures the sharpness of the operator at kd- 

To find the operator satisfying these conditions we use a multiparameter optimization of the 
coefficients mj in order to maximize D'(kd) within the constraints. We have developed a novel 
Monte Carlo routine to perform the optimization. It evolves the solution by testing randomly 
selected nearby points, selecting the best among them and iterating with a search radius sensitive 
to the speed of improvement of the solution. Because different parameters have widely varying 
ranges, we use a logarithmic sampling distribution. 

In Figure [2] we show the resulting optimized diffusivities for 5 = 4 and 5 = 5. The coefficients 
for these operators are given in Table El Comparing the 5 = 2 hyperdiffusivity to the optimized 
5 = 4 operator gives another example of the benefit of taking advantage of two free parameters. 

4. Timestep-friendly Laplacian diffusion 

Some applications require such large physical diffusivity that it becomes the dominant con- 
straint on the timestep. Examples include magnetized turbulent flows with separated viscous and 
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Table 1: Coefficients for timestep- friendly Laplacian diffusion operators 







m 




m 2 


m 3 


ffl4 


m 5 


a 2 


2 nd order 


-2 


1 










a 2 


Qth Qrc J er 


-2.72222 


1.5 


-0.15 


0.01111 






kd,-- 


=1/2 


-1.514721 


0.5692471 


0.2524535 


-0.0643401 






k d -- 


=3/8 


-0.9148711 


0.2609054 


0.2059354 


-0.0094052 






kd z 


=1/4 


-0.4334820 


0.0752718 


0.0696988 


0.0717703 






kd- 


=7/32 


-0.2733333 


0.0179722 


0.0172444 


0.1014500 






kd~- 


=3/16 


-0.2679784 


0.0393498 


0.0283651 


0.0304570 


0.0358173 




kd- 


=9/64 


-0.1317918 


0.0097316 


0.0094934 


0.0088189 


0.0081529 


0.0296992 


kd~- 


=3/32 


-0.0585600 


0.0029213 


0.0030688 


0.0029842 


0.0027554 


0.0026748 



Note. — The coefficients for a finite difference operator for d^- The 2nd and 6th order entries are for a polynomial 
fit on a radius 1 and 3 stencil, respectively. The "fed" entries are the tuned diffusion operators discussed in §[4]These 
functions are shown in Figure [T] The last entry, for ka — 3/32, additionally has mg = 0.0024189, mj = 0.0024876 
and m 8 = 0.0099690. 



Table 2: Radius-3 tuned diffusion values 



0(1/4) 


0(1/3) 


0(1) 


.124 


.328 


1.5 


.116 


.312 


2.0 


.101 


.281 


3.0 


.086 


.250 


4.0 


.055 


.187 


6.0 


.025 


.124 


8.0 



Note. — Values of the diffusion function D(k) for k = 1/4, k = 1/3, and k = 1, for the sequence of radius-3 
timestep-friendly hyperdiffusion functions with varying values of D(l). 



Table 3: Coefficients for timestep-friendly hyperdiffusion operators 

niQ mi m 2 ffl3 1714 m 5 

a 4 , 4th order 1.500000 0.250000 ~ ~ ~ 

a 4 , 6th order 1.750000 -1.218750 0.375000 -0.031250 

Tuned, stencil 4 1.231682 -0.775549 0.074534 0.126481 -0.041307 ••• 

Tuned, stencil 5 0.911455 -0.511864 -0.033699 0.094005 0.010574 -0.014743 
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resistive scales so that the magnetic Prandtl number 

Pm = ~ (18) 

V 

is far from unity, where rj is the resistivity; and turbulence with a passive scalar such as temperature 
that diffuses at a scale different from the viscous scale so that the Schmidt number 

S^- (19) 

is far from unity, where k is the diffusivity of the passive scalar. In this case, the physical diffusivity 
operators can also be adjusted so as to cause less harm to the timestep. 

The procedure that we used to generate coefficients for flat Lagrangian diffusion operators is 
to specify a value for kj, and then constrain D so as to not further increase beyond its value at the 
diffusion scale. Specifically, we set D(k) < (vrfcrf) 2 for all k. Within this constraint, we minimize 
the value of (irkd) 2 — D(k) over k > k&. For a radius S = 3 stencil, this procedure works for kd as 
low as 7/32. Any lower than that and D(k) < (irk^) 2 cannot be satisfied without D(k) taking on a 
dangerously small value for some k > k^, or even worse, becoming negative. However, increasing 
the stencil size beyond 5 = 3 allows for flat diffusion operators with successively lower values of kj. 
Such operators are shown in Table [H and Figure [1] for S = 3. 



5. Summary 

We have presented techniques for customizing diffusion filters with the goal of either decreasing 
low-A; diffusion, or maximizing the timestep, or some combination of both. We have given concrete 
examples that cover the commonly encountered cases, but since the requirements for diffusion can 
be problem dependent, we also emphasize techniques for customizing general diffusion filters. 

Turbulent flows offer a major example of the need for careful choice of the magnitude of either 
physical diffusivity or hyperdiffusivity. The relevant magnitude is that at the di ffusion scale vp(kd ) , 



where k^ is chosen to match the spectral resolution of the numerical scheme (jMaron et all 12008) . 
There it must be large enough to absorb the energy from the turbulent cascade reaching that scale. 
The value of v is generally set empirically to satisfy this requirement. 

The techniques developed here can also be applied to models with Prandtl and Schmidt num- 
bers that are large or small compared to unity, as well as models with diffusive chemistry. 



We thank J. S. Oishi for useful discussions. We acknowledge partial support of this work by 
NSF grant AST06-12724, and NASA grant NNX07AI74G. 



-10- 



REFERENCES 

Brandenburg, A. 2003, in Advances in nonlinear dynamos, eds. A. Ferriz-Mas & M. Nunez, (Taylor 
& Francis, London), 269 

Brandenburg, A. & Dobler, W. 2002, Comput. Phys. Commun. 147, 471-475 

Borue, V., & Orszag, S. A. 1995, Phys. Rev. E, 51, R856 

Canuto, C, Hussaini, M. Y., Quarteroni, A., Zhang, t. A. 1987, "Spectral Methods in Fluid 
Dynamics," Springer-Verlag 

Dedner, A., Kemm, F., Kroner, D., Munz, C.-D., Schnitzer, T. & Wesenberg, M. 2002, JCP 175, 
645 

Dunigan, T. 2004, www.csm.ornl.gov/ dunigan 

Frigo, M., & Johnson, S. G. 1998, ICASSP conference proceedings, vol. 3, pp. 1381-1384, "FFTW: 
An Adaptive Software Architecture for the FFT." 

Lele, S. K. 1992, J. Comp. Phys. 103, 16 

Londrillo, P. & Del Zanna, L. 2000, ApJ 530, 508 

Maron, J., Mac Low, M.-M., & Oishi, J. S. 2008, ApJ, 677, 520 

Maron, J. & Goldreich, P. 2001, ApJ, 554, 1175 

San Diego Supercomputing Center technical specifications, 2003, 
www. sdsc . edu/PMaC/Bencnmark/maps_ping/maps_ping_results .html 

Stone, J. & Norman, M. 1992a, ApJS, 80, 753 

Stone, J. & Norman, M. 1992b, ApJS, 80, 791 

Tamm, C. & Webb, J. 1993, JCP 197, 262 

Toth, G. 2000, JCP, 161, 605 



This preprint was prepared with the AAS IATj^X macros v5.0. 



- 11 - 



7 



4 



i 



i 



1 


1 1 1 


1 




— 

- 
- 

— 






(rrk) 2 / _ 

/ 
















/ 6 th order 








^^^^ 

x — 

2 nd order 














k d =3/8~- 












1 1 1 1 









k 



i 



Fig. 1. — Values of d^cos^nkx) = (irk) 2 , and the different finite difference operators having coeffi- 
cients listed in Table [TJ Two examples of polynomial-based difference operators with stencil radii 
S = 1 and highest order two and S = 3 and highest order six are shown, along with three examples 
of operators tuned with different choices of the beginning of the diffusive range kj. 
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Fig. 2. — Timestep-protecting diffusion functions are compared to standard hyperdiffusivities. Each 
of these functions approaches k = approximately as k , and each has been normalized so that 
.0(1/2) = 1. The function labeled "2" is the radius-2 stencil finite-difference formula for d A , or in 
other words, a V 4 hyperdiffusivity. The function labeled "3" is the radius-3 stencil finite-difference 
formula for <9 4 . Since it is higher-order, it more faithfully represents the function k 4 for large k. 
However, this is a liability for the timestep because it is more vulnerable to a diffusive timestep 
instability at k = 1 than the radius-2 function. The function labeled "4" has a different objective. 
It is a radius-4 stencil finite-difference formula for d 4 , but instead of using the extra degrees of 
freedom to represent k A at higher order, they are used to minimize D{k) for k > 1/2. The function 
labeled "5" has the same goal as "4" implemented with a radius-5 stencil. 
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Fig. 3. — The family of diffusion functions D{k) for 1 < D(l) < 8 for a stencil with radius S = 3, 
represented with solid lines. The dotted line shows the 5 = 1 operator for d 2 ; the addition of 
the two additional parameters in the S = 3 operator allowed reduction of diffusivity at both high 
and low wavenumber while maintaining the normalization at kd- The D(l) = 4 line is identical to 
the standard S = 2 hyper diffusivity. The functions are all continuous through D(kd), so the most 
diffusive above kd are the least diffusive below it. 



